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AN ALGORITHM FOR THE WEIGHTING MATRICES IN THE SAMPLED-DATA 
OPTIMAL LINEAR REGULATOR PROBLEM 

Ernest S. Armstrong and Alper K. Caglayan 
Langley Research Center 

SUMMARY 

The sampled-data optimal linear regulator problem provides a means whereby a 
control designer can use an understanding of continuous optimal regulator design to pro- 
duce a digital state variable feedback control law which satisfies continuous system per- 
formance specifications. A basic difficulty in applying the sampled-data regulator theory 
is the requirement that certain digital performance index weighting matrices, expressed 
as complicated functions of system matrices, be computed. This report presents infinite 
series representations for the weighting matrices of the time -invariant version of the 
optimal linear sampled-data regulator problem. Error bounds are given for estimating 
the effect of truncating the series expressions after a finite number of terms and a method 
is described for their computer implementation. A numerical example is given to illus- 
trate the results. 

INTRODUCTION 

Optimal linear quadratic regulator theory, currently referred to as the Linear 
Quadratic Gaussian (LQG) problem (ref. 1), has become one of the most widely accepted 
methods for determining optimal control policy. In the continuous dynamics version of 
the LQG problem, the system to be controlled is modeled as a system of first-order 
vector-matrix ordinary differential equations linear in the state and control variables. 

An optimal linear state variable feedback control law is obtained from the minimization 
of an integral performance index whose integrand is composed of weighted quadratic terms 
in the state and control. The user adjusts the weights to cause the closed loop response of 
the dynamic system to satisfy required design specifications. (See refs. 2 and 3.) A dis- 
crete dynamics analog of the continuous LQG problem also exists (ref. 4) in which the state 
equations are first-order vector-matrix finite-difference equations linear in the state and 
control variables evaluated at distinct time points. Both the continuous and discrete ver- 
sions of the LQG problem provide the control designer with a rigorous tool for developing 
linear state variable feedback control laws for multi-input multi-output dynamical systems. 

Virginia Polytechnic Institute and State University, Blacksburg, Va. 


In many control applications the physical system to be controlled is modeled by using 
continuous dynamic equations, and the control design specifications are given in terms of 
response criteria for continuous systems, but the control law is to be implemented in a 
digital fashion. A typical example is aircraft flight control applications where an onboard 
digital computer is used to generate the controller commands. In such cases a problem 
arises as to how to find a digital control law which achieves continuous design specifica- 
tions. One approach, which is in the context of the LQG methodology, is the Optimal 
Sampled-Data Regulator (OSR) problem (refs. 4 and 5). The OSR problem is a continuous 
dynamics LQG problem with one additional constraint. The operational time interval is 
divided into segments (sampling intervals) within which the control function takes on con- 
stant values. In other words, the OSR problem is a continuous LQG problem where the 
control is required to be piecewise constant. The OSR problem can easily be transformed 
into a discrete LQG problem (ref. 6) and solved by using standard LQG solution algorithms. 
The weighting matrices in the continuous OSR performance index are chosen to achieve the 
continuous design specifications as if a continuous LQG solution is to be implemented and 
then transformed into equivalent weighting matrices for use in the discrete LQG problem 
computations. The transformation equations produce extra control-state cross-product 
terms in the digital performance index and off-diagonal entries in the weighting matrices 
for the other terms which have the effect of weighting the state and control variables at 
points within the sampling interval as well as at the sampling instants. The OSR approach 
provides a means whereby the users' intuitive understanding of continuous design problems 
and LQG techniques can be applied to produce digital feedback control laws which satisfy 
continuous performance specifications. 

A primary difficulty in implementing the OSR problem is that the transformation 
equations defining the digital weighting matrices are complicated functions of the continu- 
ous system state transition matrix and other matrices appearing in the continuous formu- 
lation. In fact, the computation of these matrices for high order systems often discourages 
the use of the OSR methodology. In this report this difficulty is eliminated for the time- 
invariant version of the OSR problem by presenting a general purpose, numerically attrac- 
tive, easy-to-implement algorithm for computing the discrete weighting matrices. For the 
time -invariant OSR problem, the transformed weighting matrices are noted to be analytic 
functions of the sampling interval. Applying the analyticity property yields infinite series 
expressions for the transformation equations. Error bounds are presented that show the 
effect of truncating the series solutions after a finite number of terms and a technique is 
recommended for implementing the algorithm on a digital computer. Finally, the results 
of this paper are demonstrated with a numerical example. 
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SYMBOLS 


A 

B 

Cq(J) 

D(j) 

E(f) 

EqW, 

F 

fl, . . 
G(j) 

H(a,b) 

I 

ij 

K 

L(t) 


open-loop response matrix for continuous linear system 

control effectiveness matrix for continuous linear system 

matrix defined by equation (21) 

matrix defined by equation (59) 

matrix defined by equation (38) 

matrix defined by equation (43) 

2« + 1)||a|1(«*1)(mj)(«+2)||q ||e2«AllAti 

— where f is a natural number 

(£ + 2 ) ! 

^w('^ 2 ).E^(n 3 ^ error after truncating Q^, W, and series after nj, 

ri 2 , and ng terms, respectively 

exponential of argument c 

matrix defined by equation (B2) 

. fg functions defined by equations (B3) to (B8) 

matrix defined by equation (62) 
pa-b 

= \ e dr B where a and b are parameters 

n X n identity matrix 
indices 

positive integer satisfying condition (71) 
matrix defined by equation (58) 
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N 


n 


Qc 

Qd(i) 

%{^h) 

Rd(i) 

r 

S 

T 

t 


u(t) 

V(t) 

W(i) 


number of sampling points 
order of matrix A 

number of terms at which summation index for Qd, W, and Rd series are 
truncated, respectively 

constant nonnegative definite symmetric matrix 

matrix defined by equation (6) or (11) 

matrix defined by equation (11) 

constant symmetric positive definite matrix 

matrix defined by equation (8) or (13) 

matrix defined by equation (13) 

matrix defined by equation (54) 

number of columns in matrix B 

constant symmetric nonnegative definite matrix 

final value of time 

time variable, 0 S t S T 

ordered set of points within jjO,T^ so that tg < < . . . < tjj 

T, final time 
0, initial time 
control vector 

matrix defined by equation (34) 
matrix defined by equation (7) or (12) 
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W^Atj) matrix defined by equation (12) 

X(r) matrix defined by equation (17) 

x(t) state vector 

Y(j) matrix defined by equation (25) 

z,T,w integration variables 

A = At^/2^ 

Atj^ sample interval, = t^^j - t^^ 

II II a matrix operator norm such that ||C|| = ||C'|| for a matrix C 

Subscripts: 


c 

continuous dynamics 

d 

discrete dynamics 

i 

ith stage 



N 

final stage 



o 

initial stage 


Q 

related to. 

Qd 

computation 

R 

related to 

Rd 

computation 

R 

related to 

Rd 

computation 

W 

related to 

w 

computation 


Primes denote a matrix transpose, A dot denotes differentiation with respect to t. 
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THE TIME-BSrVARIANT SAMPLED-DATA OPTIMAL 
LINEAR REGULATOR PROBLEM 


If the linear time -invariant system 


X = Ax + Bu ^x(O) = x^^ (1) 

is given where the matrices A and B are constant and of order n x n and n x r, 
respectively, the sampled-data regulator problem occurs when the control function u(t) 

(0 = t £ T) is required to minimize the functional J 


J = x’(T) S x(T) 


+ 



x'(r) Qp 


x(r; + u’(r) R^u(t) 


dr 


( 2 ) 


subject to the restriction that u(t) be constant over subintervals tg < t-j^ < , . . < tj^ of 
[O, where 


and 


That is, 


u(t) = u^t.^ ^t. = t<L^j where i = 0, 1, . . ., N - l) (3) 

The matrices S, Q^, and R^ of equation (2) satisfy 
S - S’ 2 0 


and 


R = R' > 0 
c c 

Applying condition (3) to equation (1) gives 
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(4) 


x(t) = dT B u(tj 



Substituting equation (4) into equation (2) yields 



+ u’(ti) Rd(i)u(t.) 


with 


Q,(i) = j; ; 


tl+l A'(r-ti) A(r-t,) 


e ' -'Q^e ' "'dr 


and 


W(i) ^ r*i+l 
2 . 


R, 




dr 


H 


(‘.‘i) = i;‘ 


,A(t-T) g 


In terms of 


Ati = t.^l-ti 


(5) 


( 6 ) 

(7) 


( 8 ) 


(9) 


( 10 ) 


equations (6) to (9) become 
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( 11 ) 


«d<‘) = «d('^‘i) = c*' 

.All A-r 

2 2 Jq c 


c 


„At , 
e dr 


H(r,0) dr 


( 12 ) 



(13) 


(14) 


The sampled-data regulator problem then becomes the standard discrete linear 
quadratic regulator problem of choosing u^t^j (i = 0, 1, . . N - 1) to minimize equa- 
tion (5) subject to equations (11) to (14) and 


(^i+l) = ® '^(0 + H(At.,o) n{t.) 


i^o) = ^o 


(15) 


the solution of which can be obtained through well-known methods (ref. 4) once the matri- 
ces W^At^^, and R^^At^) are computed. The following sections of this paper 

present an algorithm for efficiently computing Q^^At^^, W ^At^ and H(j(^^i) for 3-rbi- 
trary (constant) A, B, Q^, R^, and At^. 

EVALUATION OF 
Recall from equation (11) that 



X(r) dr 


(16) 


where 


X(r) = eA’"Q^ e^" 


(X^) 
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is analytic in t and hence can be expanded in a convergent Taylor series about r = 0 


X(r) = ^ 
j=0 


djx(r) 

djj 


T=0 


j! 


Noting that 


(18) 


= A'X(t) + X(r) A 
dr 


(19) 


gives 


d^^^X(T) ^ djx(r) ^ djx(T) ^ 
dT^'''^ dr^ dr^ 

and equation (18) can be rewritten as 


(j = 0, 1, . . .) (20) 


X(r) , I 


Cq(J) 




3=0 


where 


( 21 ) 


CqO + 1) = A' Cq(J) + Cq(]) a 


From equations (16) and (21), 



(Cq(0) = (22) 


(23) 


Q, 




Cq(J) 




i+1 


i=o 


(j + 1) ! 


For computational purposes, let 


(24) 
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(25) 


Y(j) = Co(j) 




j+1 


(rri)7 


From equation (22) 


/A’ At.\ /A At.\ 

Y(i + 1) = (— -^jY(j) + ^1 


(j = 0, 1, . . Y(0) = At.j 


(26) 


whereby 


Qd(Ati) = Y(0) + Y(l) + . . . 


(27) 


In practice, equation (24) or (27) is truncated after a finite number of terms, for 
example, n^. From equation (22), 


Cq()+ l)||s 21|A|| ||Cq(1)|| 


(||Cq(0) 


Q, 


(28) 


whereby 


CqO) 


g 2^ Ia 


Q. 


(j = 0, 1, . . .) (29) 


Also, with equation (28) 

i+1 


Cq® 




j=nj+l 


(j + 1) ! 


Q. 


2^ 


' ("i) 


i+i 


j = n. + l 


(i + 1) ! 




(Equation continued on next page) 
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(" 1 + 2)1 

2("i+‘)||a1|("i"')(a,j(''i"^)||q 


i=o 


^ j + + 2 U 


(^nj + 2^! 


cil ^2|lAllAti 


(30) 


Therefore, by denoting the terms in after term n^^ by E^n^^, 

(^‘i) = 


with 


EQ(ni) < E(ni) = 




("1 ■*■ 2 ) 


The function provides an upper bound on the error obtained by truncating the 

Qd('^ti) series after n^^ terms. 

EVALUATION OF W^At.) 


Since, from equation (12), 
W^At.^ ^AL 


f 'v(t) 
*^0 


dr 


(33) 


with 


V(t) = e^^’ *^c I ^ ® 


(34) 


analytic in t, V(t) can be written as 
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V(r) = ^ 
J=0 


djy(T) 

dT^ 


II 

j! 


T=0 


( 35 ) 


Differentiating V(t) from equation (34) gives 


= A' V(t) + X(t) B 
dr 

where X(r) is given by equation (17). Recursively, 


V(T) = A- siva + aM) B 
dr^^ dri dr^ 


Defining 


CwO) = 

^ drJ 


T=0 


gives, from equation (34), 

00 

v(t) = y c (j) ^ 
^ w ]! 


i = 0 


(] = 0 , 1 , . . .) 


(36) 


(37) 


(38) 


(39) 


where for j = 0, 1, . . ., 


C^(j + 1) = A' Cw(j) + CqO) B (Cw(0) = 0; Cq(0) = 


Equation (33) yields 
W 


Ati) pAt^ 

= \ V(T) dT 
2 Jq 


(40) 


(41) 
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or 



I c„(J) 


j=l 



(j + 1) ! 


For computational purposes, let 


D(j) 


Cw(i) 


kT 

(j + 1) ! 


(42) 


(43) 


Then for j = 0, 1, . . 


/A’ At.\ /B At. 

D(i + 1) = + Y(j)‘ 1 


(^) {''») = »! V(O).Q^Atj) 


(44) 


where Y(j) is computed from equation (26). In terms of D, 
W^At.) = 2 |p(l) + D(2) + .. 2 


(45) 


The computation of an upper bound on the error in truncating equation (42) after a 
finite number of terms, for example, j = n 2 § 1, is somewhat more complicated than the 
corresponding computation for (At.) but can still be carried out by using the same 
general approach. From equations (29) and (40), 


Cw(j + 1) 


= A 


Cq(J) 




Cw(i)| 




Rearranging equation (46) and assuming that j|A|| o yields 



-j 




-j+1 






A 


Cw(i + 1) 

- 

A 


C-^(j) 

5 2^ 

Q 

^c 


B 


Summing inequality (47) from j = 0 to j - 1 yields 


A 


-j+1 


C^(j) 



Cw(0) 





(46) 


(47) 
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or 




. llQjl ||b|| ||A|p-*(2i - l) 


Writing 


W 


(*‘i) 


= 2 


Cwd + dK)'"'^ ^ ^ Cwd 


2 f^o (i + 2 ) ! 


(j + 2)! 

]=n2+l 


(48) 


(49) 


and employing equation (48) gives 

CwO + 


s 

j=no+l 


(j + 2)! 


j=n2 + l 


||Cw(j + 1)|| ('^^i)^ 

(j + 2) ! 


Q, 




z 

j=n2 + l 


(j + 2)! 


i xio + 2 ^ 

fn.,+2')/ .(n^ + s) 

2^ 2 Vll 



(„2 3) ! 





21|a|I At. 


IIbH 

INI 


(50) 


Thus 


W At. 


^2 ■ - ' 


ii z J. / \ 

- “ 2 (I H- i Ti — 

i=i 


(51) 


where 


EwW =®('^2+l)|^ 


l|A|l 


(52) 


for A #0 and n^ S 1. 


EVALUATION OF (^t- 


The relation R^(^AtJ can be written as 


14 


( 53 ) 


where 

(^^i) = Jo 


Define L(t) as 



(54) 


(55) 


from Which 

= B'V(r) + V'(t) B 

with V(t) defined from equation (34). Recursively, 

L(t) = B’ ^ (t) + (t) B (i = 0, 1, . . .) 

dri+1 drl drl 

The function L(r) may alternately be expressed as 


ji 

i! 


With 




djL(T) 


1 = 0 


drJ 


Cp(j) 


d^L(r) 

dT^ 


T=0 


(56) 


(57) 


(58) 


(59) 
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equation (57) gives 


C^ii + 1) = B'C^(j) + C^(j) B 


(60) 


with 


Cp(0) = Cp(l) = 0 


Then 


= i L(r) dr = ^ Cj,(j) 
^ j=2 


(61) 


For computational purposes, let 

G(i) = y — 

(i + 1) ! 


(62) 


with 


B' At.\ B At.' 

G(j + 1) = I i D(j) + D'(j) i 

• - \j + 2y 


\i + 2, 


(i = 0, 1, . . .) 


(63) 


from which 


Rd(Ati) = G(2) + G(3) + . 


(64) 


By following the same procedure as in the previous sections, an upper bound on the 
error in truncating equation (61) after j = n^ ^ 2 terms can also be found as follows: 



j=0 


Cr(J + 2) 


(j + 3) ! 



(65) 
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( 66 ) 




and 



(67) 


SUMMARY OF RESULTS 


In the previous sections it has been established, for ||A|| i= 0, that 




e^^ dT = Y(0) + Y(l) + . . . + Y(n^) 



Y(j + 1) = 


fA' At, 


3 + 2 


i Y(j) + Y(j) 


At. ^ 

j + 2y 


(y(0) = 






2} I a) i At. 
e 


W 


^At.^ pAtj^ *, pT . 

^2 ~ Jo ® ^ Jq- ® B dz dr = D(l) + D(2) + . . . + D(n 2 ) + Ew(n 2 ) 


D(j + 1) = ^ D(j) + Y(j) (^(0) = 0; Y(0) = At. 


Ew{n2)| <E(n2 


+ 1 


liAll 
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= G(2) + G(3) + . . . + G(ng) + E^(n3) 



The absolute error bounds for EQ^n^^, E^^n 2 ^, and E^^ng^ can, in each case, be 
written in terms of the error function E. Relative accuracy requirements, however, 
generally require different values for nj, n 2 , and n 3 « 

For the trivial case in which ||A|| = 0, 


!d("l) 

= Qc 

(68) 

/(At,) 

= QcB(At.)' 

(69) 


(al)^ 


d(«i) 


(70) 


COMPUTER IMPLEMENTATION 


The infinite series representations summarized in the previous section, when 
numerically evaluated, yield Qcj(At^), W^At^^,and R^^At-^ for general A, B, Q^, 
and R^ matrices and sampling intervals At^^. In practice, numerical difficulties can 
occur unless the matrix A is preconditioned to avoid large numbers of terms in the 
series solutions and excessive roundoff error. One preconditioning approach and its 
implementation are presented in this section. 
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If the matrix A and sample time At^ are given, find an integer K g 0 such that 



By defining 


A = 



(71) 


(72) 


inequality (71) gives 


AA 


< 1 


(73) 


whereby the eigenvalues of AA are within the unit circle in the complex plane and 

lim (AA)^ = 0 (74) 

i— oo 

Next, replace At^ by A and compute the corresponding Q^(A), W(A), and Rjj(A) 

matrices by the series solutions. Convergence should be rapid because of conditions (73) 


and (74). 
W(A), and 
tions (A5), 

If K > 0, the required solutions for Atj^ may be constructed from 
R^(A) by applying the equations derived in appendix A as follows. 
(A7), (A8), and (A9) replace A by 2^A to give 

Qd(A), 

In equa- 


Q^(2j+^A) = Q^(2jA) + Q^{2^a) 

(75) 


w(2'’'^^a) = (i + w(2jA) + 2e^’^^"^ Q^(2jA)H(2jA,o) 

(76) 


R^(2^+^a) - 2R^(2jA) + W’(2U) H^2jA,o) + H' ( 2 ^, 0 ) 



+ H'(2jA,o)Q^(2U)H(2jA,o) 

(77) 


h(2^^^A,o) = (i + e^^^^)H(2^A,o) 

(78) 


eA2^'^^A ^ e^^^A 

(79) 
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Recursively evaluating equations (75) to (79) from j = 0 to j = K - 1 yields Q^^At^j, 

W^At^^, and R^^At^^ at the final stage. Algorithms for computing e"^^ and H(A,0) 
to initialize the recursive process are given by Kallstrom. (See ref. 7.) 

Software for the theory presented in this paper is available in the ORACLS program. 
(See refs. 8 and 9.) 


AN EXAMPLE COMPUTATION 


In this section a particular set of ^A, B, Atj^ is chosen and results from the 
algorithms presented in the foregoing sections are illustrated. The computation was 
performed by using the subroutine SAMPL of the ORACLS program (ref. 9) on a CDC 
6600 digital computer in single precision. The SAMPL subroutine employs the method of 
computer implementation described in the preceding section. Numerically, convergence 
was assumed to have occurred in the Q^, W, and R^ series when the improvement in 
the element of largest magnitude (measured relatively if the magnitude was less than unity, 
and absolutely otherwise) of each of the matrices was past the eighth significant digit. Let 



1 1 
2 1 
1 2 


and 


B = 


0 

1 

0 



0 0 
1 0 

0 1 



0 

1 


At. = 1/2 
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Taking 


m 

r II 


= max ; mc-jj 
i=l 


(1 ^ j g n) 


for some given m x n matrix C with elements c^j gives 


and 


A = A1 = 4 


B = B’ = 1 


The condition 
At 


||A|| 

gives K = 2 and A = 1/8. The results of the series computations are 


Qd(l/8) 


f (1/8) = 


0.166270 

0.0242575 

0.0242575 

0.0242575 

0.166270 

0.0242575 

0.0242575 

0.0242575 

0.166270_ 

0.0102932 

0.00142899 


0.00142899 

0.0102932 I 


0.00142899 

0.00142899J 


r 


“1 


Rd(l/8) = 


0.000798547 

0.0000827449 


0.0000827449 

0.000798547 


with truncation error bounds 

i|EQ(ni)|| < E(nj) = 7.09 X lO'^O (n^ = lo) 
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("2 = 9 ) 

("3 = s ) 


|Ewh)|l ^ + l) |{f}{ = 





10“^^ 


Relative accuracies are 




||Qd(i/8: 

'1 

||V2)| 


f(V8) 


®r("3) 



|Rd(l/8) 


< 3.30 X 10“® 
1.35 X 10“® 


< 5.03 X 10"^ 


Since K = 2, two passes through equations (75) to (77) were required which generated the 
sequence 


Qd(l/4) = 



0.482451 

0.158090 

0.158090 

0.158090 

0.482451 

0.158090 

0.158090 

0.158090 

0.482451 

0.0576453 

0.0173101 


0.0173101 

0.0576453 


0.0173101 

0.0173101 



Rd(l/4) 


0.00815427 

0.00184446 


0.00184446 

0.00815427 


Qd(l/2) = 


2.80602 

1.94688 

1.94688 


1.94688 

2.80602 

1.94688 


1.94688 

1.94688 

2.80602 
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f (1/2) = 


0.565488 

0.355069 

0.355069 


0.355069 

0.565488 

0.355069 




0.124575 

0.0628764 


0.0628764 

0.124575 


Rd(l/2) 


0.624575 

0.0628764 


0.0628764 

0.624575 


These final results agree with the closed form solutions shown in appendix B. 

CONCLUDING REMARKS 

A computational procedure for generating the weighting matrices needed in the time- 
invariant optimal linear sampled-data regulator problem has been presented. This pro- 
cedure makes use of an analytical property of the defining equations to produce general 
purpose numerically attractive infinite series expansions which can be easily summed on 
a digital computer. Error bounds for truncating the series expansions after a finite num- 
ber of terms were derived and gave good agreement with intended accuracy for the numeri- 
cal example considered. It is felt that the results of this paper eliminate a basic and 
major difficulty in optimal sampled-data regulator methodology and open the way to a 
wider application of the theory. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
November 11, 1976 
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APPENDIX A 


EVALUATION OF SAMPLED-DATA WEIGHTING MATRICES FOR 
DOUBLED VALUES OF SAMPLING INTERVAL 


For the purposes of this appendix, let the sampling interval be denoted by A and 
the weighting matrices defined by equations (11), (12), and (54) be denoted by 



(Al) 


W(A) _ 
2 



H(r,0) dr 


(A2) 



H'(t, 0) H(r,0) dr 


(A3) 


with, as before. 


H(r,0) = r 

- 'n 


e^^ dz B 


(A4) 


Here equations for Q^(2A), W(2A), and R^(2A) are derived in terms of Qjj(A), W(A), 
R^(A), e'^^, and H(A,0). These equations are used in the body of the paper in the section 
dealing with a computer implementation of the algorithms for Q,( At- ), w(At.V and 

For Q^, 


R 


Q^(2A) = e^^dr = Q^(A) + d^ 


A' A 


= Q^(A) + e" “ Q^(A) e 


AA 


(A5) 


Additionally, 
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H(t + A,0) 




T + A 


B dz = H(A,0) 


I" 

‘Jo 


A(z+A) 


B dz 


= H(A,0) + H(t,0) 


from which 


H(t+A,0) = H(A,0) + H(t,0) = H(t,0) + H(A,0) 


(A6) 


and 


H(2A,0) = ( 
For W, 


I + e 


AAj 


H(A,0) 


(A7) 


W(2A) ^ gA’r 

2 Jo 


H(r,0) dr = 


+ eA’(r+A) h( 7+A,0) di 


whereby 


^ W(A) 
2 


+ e 


A' A 




H(t, 0) dr + 



e^'^ dr H(A,0) 


W(2A) = (l + e^'‘^)w(A) + 2e^’^ Q^j(A) H(A,0) 


(A8) 


Finally for R^j, 


^ P 2 A ^ P A 

R^(2A) = H’(t,0) H(t, 0) dT = R^(A) + H’(r+A,0) H(t+A, 0) dr 


= R^(A) + H'(t,0) H(t,0) dT + e^^ dr H(A, 


0 ) 


(Equation continued on next page) 
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+ H*(A,0) \ e"^’^Q H(t, 0) dr + H’(A,0) C ^ Qc e'^^drH(A,0) 

Jq *^0 


2R (A) + h(A,0) + H’(A,0) + H’(A,0) Q^{A) H(A,0) (A9) 

Cl 2 " 


APPENDIX B 

A CLOSED FORM SOLUTION FOR COMPUTATIONAL EXAMPLE 

A closed form solution is given for the numerical example presented in the body of 
this paper. These results are useful in determining the validity of the numerical 
computations. 

For 



e^* = eh + |■(e‘^* - e*)F 

where I is a 3 x 3 identity matrix and 
F = A - I 

With equation (Bl), 



(Bl) 


(B2) 




equations (10), (11), and (12) yield 
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where 


Q„{Ati) = 


^1 + h 


^1 ^2 


^1 ^ ^21 


W(At.) = 


^3 + ^4 


f3-^^4 






% ^6 



(B3) 

(B4) 

(B5) 

(B6) 

(B7) 
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